/* This function takes two cmatrices and multiplies it. */
# include <stdio.h>
# include <stdlib.h>
# include <math.h>

/* Number of rows and columns of the left and right matrices to be multiplied. */
# define NROWSL 2
# define NCOLSL 3

# define NROWSR NCOLSL
# define NCOLSR 4

struct cmtrx
   {
    unsigned short nrows, ncols, padpti, padpte;
    short *istart, *iend; // If every element in the row is zero, istart[row] = -1, iend[row] = -2;
    double **mtrx;
   };

extern short multiply_cmtrx(struct cmtrx *mleft, struct cmtrx *mright, struct cmtrx *mres);

int main()
 {
  short i, k;
  struct cmtrx mleft, mright, mres;

  mleft.nrows = NROWSL;
  mleft.ncols = NCOLSL;
  mleft.istart = (short *)malloc(sizeof(short)*mleft.nrows);
  mleft.iend = (short *)malloc(sizeof(short)*mleft.nrows);
  mleft.mtrx = (double **)malloc(sizeof(double *)*mleft.nrows);
  for(i = 0; i < mleft.nrows; i++) mleft.mtrx[i] = (double *)malloc(sizeof(double)*mleft.ncols);

  mright.nrows = NROWSR;
  mright.ncols = NCOLSR;
  mright.mtrx = (double **)malloc(sizeof(double *)*mright.nrows);
  mright.istart = (short *)malloc(sizeof(short)*mright.nrows);
  mright.iend = (short *)malloc(sizeof(short)*mright.nrows);
  for(i = 0; i < mright.nrows; i++) mright.mtrx[i] = (double *)malloc(sizeof(double)*mright.ncols);

  mres.nrows = mleft.nrows;
  mres.ncols = mright.ncols;
  mres.mtrx = (double **)malloc(sizeof(double *)*mres.nrows);
  mres.istart = (short *)malloc(sizeof(short)*mres.nrows);
  mres.iend = (short *)malloc(sizeof(short)*mres.nrows);
  for(i = 0; i < mres.nrows; i++) mres.mtrx[i] = (double *)malloc(sizeof(double)*mres.ncols);

     mleft.istart[0] = 1;
     mleft.iend[0] = 2;

     mleft.istart[1] = 0;
     mleft.iend[1] = 0;

     mright.istart[0] = 2;
     mright.iend[0] = 2;

     mright.istart[1] = 0;
     mright.iend[1] = 3;

     mright.istart[2] = 0;
     mright.iend[2] = 1;


     mleft.mtrx[0][0] = 1;
     mleft.mtrx[0][1] = 2;
     mleft.mtrx[0][2] = -sqrt(5.0);
     mleft.mtrx[1][0] = 1;
     mleft.mtrx[1][1] = -sqrt(5.0);
     mleft.mtrx[1][2] = -sqrt(5.0);

     mright.mtrx[0][0] = 1;
     mright.mtrx[0][1] = -sqrt(5.0);
     mright.mtrx[0][2] = -sqrt(5.0);
     mright.mtrx[0][3] = -sqrt(5.0);

     mright.mtrx[1][0] = 1;
     mright.mtrx[1][1] = 0;
     mright.mtrx[1][2] = 0;
     mright.mtrx[1][3] = 3;

     mright.mtrx[2][0] = 1;
     mright.mtrx[2][1] = 2;
     mright.mtrx[2][2] = -sqrt(5.0);
     mright.mtrx[2][3] = -sqrt(5.0);

  printf("\nLeft Matrix:\n");
  for(i = 0; i < mleft.nrows; i++)
	{
	 for(k = 0; k < mleft.ncols; k++) printf("%E\t", mleft.mtrx[i][k]);
	 printf("\n");
	}

  printf("\nRight Matrix:\n");
  for(i = 0; i < mright.nrows; i++)
	{
	 for(k = 0; k < mright.nrows; k++) printf("%E\t", mright.mtrx[i][k]);
	 printf("\n");
	}

  if( multiply_cmtrx(&mleft, &mright, &mres) < 0) return -1;

  printf("\nMultiplied Matrix:\n");
  for(i = 0; i < mres.nrows; i++)
	{
	 for(k = mres.istart[i]; k <= mres.iend[i]; k++) printf("%E\t", mres.mtrx[i][k-mres.istart[i]]);
	 printf("Range = %d~%d\n", mres.istart[i], mres.iend[i]);
	}

/* Freeing the allocated memory. */
  for(i = 0; i < mres.nrows; i++) free(mres.mtrx[i]);
  free(mres.mtrx);
  free(mres.istart);
  free(mres.iend);

  for(i = 0; i < mright.nrows; i++) free(mright.mtrx[i]);
  free(mright.mtrx);
  free(mright.istart);
  free(mright.iend);

  for(i = 0; i < mleft.nrows; i++) free(mleft.mtrx[i]);
  free(mleft.mtrx);
  free(mleft.istart);
  free(mleft.iend);

  return 1;
 }
